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NUMERICAL EXPERIMENTS WITH FLOWS OF ELONGATED GRANULES 
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14 Cromwell Court 
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and 

D.E. Brewe 

Propulsion Directorate 
U.S. Army Aviation Systems Command 
Lewis Research Center 
Cleveland, Ohio 44135 

Theory and numerical results are given for a program simulating two-dimensional granular flows 
(1) between two infinite, counter-moving, parallel, roughened walls, and (2) for an inf initely-wide 
slider. Each granule is simulated by a central repulsive force field ratcheted with force restitution 
factor to introduce, dissipation . Transmission of angular momentum between particles occurs via Coulomb 
friction. The effect of granule hardness is explored. Gaps from 7 to 28 particle diameters are 
investigated, with solid fractions ranging from 0.2 to 0.9. Among features observed are: slip flow at 
boundaries, coagulation at high densities, and gross fluctuations in surface stress. 

A videotape has been prepared to demonstrate the foregoing effects. 


INTRODUCTION: 

The purpose of this work is to show how 
computer simulations of granular flow can 
assist in the understanding of third-body load 
support from solid particles. In a literature 
review conducted several years ago [1] we 
found that conditions of tribological interest 
most frequently begin with particle densities 
greater than those encompassed by the then- 
current extensions of classical kinetic theo- 
ry. It was concluded that "simulation programs 
are a must for third-body tribology." This 
paper reports on efforts to follow our own 
recommendation. 

In fact, since our review paper, comput- 
er simulations have come into increasing use 
[See ref. 2 tor a recent extensive review]. 
Much of the work has been directed to the 
mutual confirmation of analytical and numeri- 
cal methods. This confirmation has been sought 
through the exploration of regions of overlap- 
ping applicability between simulations and 
earlier analytical results. In particular, 
constitutive equations have been examined for 
moderate-density flows of infinite extent, 
with most of the numerical experiments being 
performed in three dimensions, and with sub- 
stantial computer-times being registered. 

If suitable constitutive equations can 
be derived, then appropriate boundary condi- 
tions must still be found. Work toward deter- 
mination of the latter has not been nearly as 
extensive. In addition, for our purposes, 
there are at least two other obstacles to 
overcome. First, in the case of narrow gaps, 
it is possible that no internal region exists 


which can be suitably approximated by the 
constitutive equations for an infinite medium. 
Second, there must be means for prediction of 
the large stress and velocity fluctuations 
which occur. Our purpose in mentioning these 
difficulties is simply to explain our own 
choice for a first area of investigation. 
Certainly for those applications where it is 
applicable, a successful continuum approach 
should yield a substantial reduction in the 
time required for performance calculations, 
and a greater understanding of the associated 
flow phenomena . We do recognize the importance 
of this approach, and hope subsequently to 
participate in its implementation. 

This paper describes our efforts to 
develop simulation programs for two-dimension- 
al flows. The adoption of this 2-D assumption 
is tantamount to the assumption of an elongat- 
ed particle shape , with axis transverse to the 
flow direction. At this time it is unclear 
whether typical particles are better approxi- 
mated by cylinders or by spheres. In the case 
of wear particles, there is some support for 
the first alternative. For us, the lower 
computer time provided by the 2-D assumption 
was decisive. A much wider range of operating 
parameters can thereby be explored, and compu- 
tations are feasible even on a desktop comput- 
er. Figure 1 provides a "snapshot" of one flow 
configuration extensively examined. It con- 
sists of two parallel, counter-moving walls of 
infinite lateral extent separated by a granu- 
lar medium. In the case of a constant-property 
Newtonian fluid, the result would be a Couette 
flow. 
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Figure 1 Granular Flow Field Between Infinite 
Counter-Moving Roughened Parallel Walls 
90% solids 



Figure 2 Geometry of Particle-Particle 
Interaction 


I = Ma 2 /2 


Through personal communication we have 
recently become aware of some research very 
similar to our own, conducted simultaneously, 
but independently of us, by S. B. Savage and 
R. Dai [3,4] (hereinafter referred to as S&D) . 
They have treated bounded simple-shear flows 
of smooth, inelastic spheres. Their "walls" 
consisted of massive spheres of the same 
diameter, spaced to mimic the interior parti- 
cle distribution. Two different models for 
interaction were treated, with similar results 
obtained from each. In cne model, the spheres 
were treated as rigid, with a coefficient of 
velocity restitution used to introduce dissi- 
pation. In the other model, "soft" spheres 
were used, with dissipation introduced via the 
Walton "ratchet" model (as also used in this 
paper). We shall report on some of their 
findings at the time we discuss our own re- 
sults . 

At the outset, we wish to advise the 
reader that we do not expect that the kind of 
investigation reported here will lead to 
numerically correct predictions of real-life 
physical systems. Rather, we anticipate that 
the results will be semi-quantitative, capable 
of providing important clues to behavior, and 
a foundation for later empirical representa- 
tions of real-life experiments. 

THEORY 

The analysis here will be developed for 
dissimilar cylinders in contact, although all 
the reported computations will be for granules 
of uniform size. 

Geometry of Collision : 

A granule is taken to be an elongated 
cylinder of length, L, radius, a, and density, 
Pp . The mass of such a cylinder is: 

M = p P n a 2 L 

And the moment of inertia of such a granule 
is: 


Figure 2 shows the interaction of two 
granules, I and J, with radii of ai and aj , 
respectively. For simplicity, we adopt the 
notation : 

= T) - n = et 6xi + §2 5x2 

r = I o£| = separation of centers 

n t = 5r/r = unit normal to point 
of contact 

Expression for Normal For ce of Interaction : 

The normal separating force exerted by 
granule J upon granule I (each of length, L) 
is taken to be: 

Fij,n = — ni G (am / aw) $ ( r / am ) fac 

Here "G" is a constant force magnitude charac- 
teristic of the granular material, "aw" is the 
radius of the wall granules, "am" is the mean 
radius of the two interacting granules, $ is 
a function to be discussed in a later section, 
and fac is a factor accounting for energy loss 
during collision. 

During the separation following a colli- 
sion, fac reduces the repulsive force between 
the granules. Separation is taking place 
whenever : 

5V> 5r > 0 

Therefore, we take fac to be: 

fac = [1 + res - (1 - res) *SGN (5V* 5r.) ] /2 

Note that fac = 1 during granule-to-granule 
approach, and fac = res during separation. 

fac is not the coefficient of velocity 
restitution , resv . The relationship is (ap- 
proximately) : 

resv = /fac 
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Expression for the Tangential Forces: 

The tangential force exerted by one 
cylinder on another during collision can only 
be represented approximately. We do not ac- 
count for the deformation which takes place at 
the surface of the cylinders. Referring to 
Fig. 2, we take the origin of the vectors nt 
and U as the point of contact. The distance 
from the origin of cylinder, I, to this point 
is: 

armij - r /2 + (a»/r)(ai - aj ) 

where: at and aj are the radii of the cylin- 
ders I and J, respectively, "r" is the separa- 
tion of their centers, and 

= ( 1 / 2 ) (aj +aj ) 

Since this expression can lead to preposterous 
numbers when the cylinder centers are very 
close, we accept values from the above expres- 
sion only if they exceed: 

armij (minimum) = ai - (aj) 2 /( 2 ai) 

Now if the cylinder I is rotating coun- 
terclockwise with angular velocity, m , then 
the local velocity at the so-defined contact 
point is: 

Ui = Vi + <ji (armij /2) ti 
The tangent vector, jti , is given by: 

1 1 = e 3 x ni - (£2 5xi - ei 5x 2 ) / r 

The tangential component of velocity of 
particle I relative to particle J at the point 
of contact is (Ui - Uj ) ■ t.i , and is 

TMT = -5V-ti + ui armij + uj armj i 

If the above component is positive, then the 
particle J tends to reduce the angular veloci- 
ty of particle I. The tangential force exerted 
by particle J upon particle I is taken to be 
proportional to the normal repulsive force, 
with a coefficient of friction, Cf . Thus: 

F, m = - SGN(TMT) ti C f G (a/a«) ${r/am) fac 

Equations of Motion: 

The equations of motion can now be 
written: 

Mi (dVi /dt) = (£ J|# + J, it ) 

Ii (dui /dt) = armij Fi M 

The summations in these equations account for 
multi-body collisions. 

Non-Dimensional Form: 

To reduce as much as possible the number 
of parameters that must be explored numerical- 
ly, we put the equations into dimensionless 
form. We adopt the following 


Unit of Distance: a* 

Unit of Time: a* /U* 

Unit of Mass: M* 

Certain irremovable dimensionless param- 
eters remain in the problem. These are: 

1. res, the coefficient of force res- 
titution 

2. Cf , the coefficient of friction 
between particles 

3. f c = Gaw / (Mw Uw 2 ) 

4. (plate separation) /2a* = W/2a* 

5. 

(separation of roughness particles) /2a* 

6. solids fraction 

7 parameters used in the dimension- 

less repulsion function, <fr(r/a*) 

Hereafter in this section, dimensionless vari- 
ables will be employed in the equations. Thus 
the new V will be V/Uw , the new u will be 
oa/Uv , etc. 

The equation for translational motion 
now becomes: 

dVi /dt = -£i Fn [m + ti SGN(TMT)Cf ] 
and that for rotation: 

dui /dt = -Ei Fri (Mw/Mi ) C/ SGN (TMT) armt j 
where: Fri = fc (Mw /Hi ) <a« /a* ) 4> (r/a« ) f ac 
and: 

fac = [1 + res - (1 -res) SGN(TMN) ] /2 
Formula for Repulsive Force: 

We revert here to dimensional symbols. 
The normal force on granule I by granule J has 
previously been written as: 

G (a* /aw )& (r/a» ) fac 

In the present study we adopt for $ a rather 
arbitrary positive function of the relative 
separation, r/ao , with the thought that any 
generalities concerning granular flow must not 
depend critically on the details of such a 
formulation. Certainly this supposition is 
true in connection with classical kinetic 
theory of gases. The function chosen for the 
present investigation is: 

0(r/a«) = exp(-ar/a«) 

We define the "force ratio", fr, as the 
ratio of the repulsive force between wall 
granules when the center-to-center separation 
is a w to that when the separation is 2a w (that 
is, when the granules are on the verge of no 
contact). Thus: 

a = log. (fr) 

This parameter, a, determines how rapidly the 
repulsive force increases with contact pene- 
tration. 
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Next we choose as a measure of the 
hardness of the granules the amount of energy 
required to force an approach from initial 
contact at 2a„ to the separation of aw. Thus: 


2a„ 

E - J 9 (r/a v ) dr 


[Ga„(/r-l) ] 
t (fr) 2 ln(fr)] 


To make this energy non-dimensional , we divide 
by the kinetic energy of a wall granule moving 
at velocity, U* , and form the "hardness ra- 
tio", hr. Thus: 

hr = 2ga » (-f-r-l) 

M V (U V ) 2 <fr) 2 ln(fr) 


In the case of two counter-moving walls, we 
take 

Uw = U(upper) - U(lower) 

As a result of the above definition: 

f = ga y _ hi (fr) 2 ln(/r) 
c M V (U V ) 2 2 (fr-1) 


The collision time may be estimated to 
be roughly twice the time to bring the parti- 
cle I to a halt. Figure 3 shows the relation 
between this collision time and the hardness 
and force ratios. When the force ratio is 
unity, this time varies inversely with the 
hardness ratio. The logarithmic scale employed 
on the chart shows that this relationship is 
also substantially correct for higher force 
ratios . 



Collision Time : 

To estimate the collision time for a 
binary encounter, let us assume that granule 
J is fixed, that granule I is moving solely in 
the x-direction with initial velocity, Uv, at 
the commencement of contact (where x is taken 
as zero) . For simplicity, we take both gran- 
ules to have the same radius, aw. The motion 
of granule I is then governed by: 

du/dt = -{Gaw/MwUw 2 } exp [- ql ( 2-x) ] 

Multiplying top and bottom of the foregoing 
derivative by "u" and integrating, we get: 

u 2 = 1 - 3 [exp(ax) - 1} 

where : 

P = (2Gaw/MwUw 3 )exp(-2a) 

The velocity of particle I is reduced to zero 
at the point: 

x #u = [In (1 + l/p)] /a 

Now with dt = dx/u, we have for the elapsed 
time after initial contact: 


{ [1 +P~Pexp (ax) ] 1/2 

This last integration can be accomplished in 
closed form. Thus: 

t = (x - 2A/a) / (1+p) 1,1 

where: 

A = ln[(l+3) t/l + { 1+3 - 3exp (ax) ) l/I ] - 
In [ (1+p) 1,1 + 1] 


Figure 3 Collision Time as Function of 
Hardness and Force Ratios 


It is important in the numerical cal- 
culations that the time-step, At, be small 
enough relative to the collision time for the 
interaction process to be adequately sampled. 
Figure 4 shows the effect of the magnitude of 
At on the shear stresses computed for a shear 
flow between two counter-moving parallel walls 
separated by a gap of 7 granules diameters 
with a 50% solids content. 

As a result of this At-test, the fol- 
lowing schedule of At was observed in all 
computations. 

HR M 

10 0.02 

5 0.04 

3 0.08 

Wall Stresses : 

The computations are necessarily per- 
formed for a finite computational volume. In 
the case of parallel walls in relative motion, 
this volume has width, W (y-direction) , length 
D (x-direction) , _and depth, L (z-direction, 
eventually taken as infinite) . Particles are 
affixed to the walls, with spacing close 
enough so that these particles absorb all the 
forces exerted on the walls. 

Let EF denote the sum of forces applied 
to the particles stretched along one of the 
walls over length, D. Then the wall stress is 


4 




u = £F/(LD) — ZM*a/(LD) 

To make the stress non-dimensional, we divide 
by pp Uw * . Thus: 

(j/ ( Pp Uv 2 ) = Zna* 2 La/ (LDU* 2 ) = 

tt [Z^a* ZOw 2 ] / {D/av ) 

The righthand side of this last equation 
involves the non-dimensional acceleration, 
aaw/tk 2 , which the wall particles would pos- 
sess if unrestrained, as well as the non- 
dimensional length, D/a w , of the wall on which 
the stress is computed. 



equisized slices or the gap. rrom rnese aata, 
densities and velocities were computed. 

The selection of a desirable time step. 
At, has been discussed. However, the total 
time to reach., a steady-state also must be 
determined. Inasmuch as the flow never attains 
a completely steady character, some judgment 
has to be exercised as to when to terminate 
the calculations. Although a more rigorous 
criterion might have been set up (with greater 
computing time required) , the decision for 
each case was made by examining the transient 



TIME 


Figure 4 Computed Results for Shear Stress 
as Function of At Employed 

NUMERICAL EXPERIMENTS 

SHEAR FLOW: 

The first numerical experiments were 
performed on shear flows between counter- 
moving parallel walls, as shown in Figure 1. 
The walls were roughened by the attachment of 
half-cylinders with a center-to-center spacing 
(unless otherwise stated) of 1.5 diameters. 
For this configuration, the gap is defined as 
the center-to-center distance between opposing 
wall granules. An infinite lengthwise dimen- 
sion was mimicked by imposition of periodicity 
in the distribution of granules. The flow 
computation was commenced by placing granules 
at random, with zero velocity, within the 
computational region. In order to attain 
equilibrium, these granules would repel one 
another, hitting the walls, and commencing the 
flow. 

Flow quantities were averaged in two 
stages. In the first stage, quantities were 
summed and averaged for 50 timesteps. The 
resultant averages were treated as macroscopic 
properties, local in time. In the second step, 
these quantities were averaged over lengthy 
time periods during which quasi-equilibrium 
obtained. The time-averaged opposing stresses 
for the second step were found to match within 
a few percent. This unforced agreement attests 
to an accurate implementation of Newton's Laws 
within the medium. The particle numbers and 
momenta were also summed and averaged for five 


Figure 5 Development of Velocities at 
1/10 and 3/10 Across Gap 

development. Figure 5 shows such development 
for the velocities at 0.1 and 0.3 across the 
gap for the flow shown in Fig. 1. In this 
particular instance, it would appear sensible 
to commence averaging for a "final” results at 
t = 200. Many particles are involved in this 
computation (1650) and computing time can 
become excessive. For example, the total 
runtime for this case was 16.5 hrs. on an IBM 
(clone) 486 33. This computer runs approxi- 
mately 36 times as fast as the original 
IBM AT. 

With so many parameters assignable for 
the investigation, some choices had to be 
made. The effect of granule hardness was 
explored with a gap of 7 diams. Figure 6 shows 
that shear stresses at densities below 60% are 
virtually independent of HR. This result is in 
keeping with classical kinetic theory of 
gases, where different molecular models tend 
to give similar results. It is, moreover, in 
agreement with the finding of S&D of similar 
results for rigid and soft granules. Above 60% 
solids, our Values rise rapidly, and at rigid- 
body complete filling (90.7%), they are pro- 
portional to HR. 

In all further discussion, unless spe- 
cifically stated otherwise, the following 
values were assigned for the numerical compu- 
tations. 

HR = 5, FR = 3, res = .8, CF = .3 
Ctr.-to-Ctr. Spacing of Wall Particles, 1.5 
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Figure 6 Shear Stress Dependence on 
Hardness Ratio and Density 



Figure 7 Velocity Distribs. for 7 Diam. Gap 
HR = 5; FR = 3; CF = 0.3; res = 0.8 



Figure 8 Particle-Count Density Distribution 
for 7 Diam. Gap 

Figures 7 and 8 show the velocity and 
particle distributions obtained for a gap of 
7 diams. Antisymmetry was forced on the veloc- 
ity profiles, with the supporting points shown 
as data. At low density, the velocity distri- 
bution is typical Couette Flow, with slip at 
the walls. With increasing density, the veloc- 
ity gradient tends to be steeper in the inte- 
rior than it is near the walls, and the slip 


approaches zero (as found also by S&D). The 
density distributions are only for the "ac- 
tive" or unattached granules. Near the wall 
the available space for such active granules 
is restricted, and the effective density is 
augmented by the inclusion of the half-cylin- 
ders (in this case making an addition of 5). 

Figures 9 and 10 show the velocities and 
densities for a gap of 28 diams. At low densi- 
ties, the Couette Flow with slip is again ob- 
served. However, in contrast with the small 
gap case, the interior velocity gradient 
diminishes with increasing granule density. In 
the high density limit, a plug is formed 
within the interior, the particle aggregate 
becoming effectively a solid. This phenomenon 
is clearly visible in Figure 1, where the 
granules are lined up within the interior, and 
ruffled in the regions adjacent to the walls. 
The density distribution shows a slight trough 
in the center of the gap. ~ 



POSITION IN GAP 


Figure 9 Velocity Dist. for 28 Diam. Gap 
HR = 5; FR = 3; CF = 0.3; res = 0.8 



Figure 10 Particle-Count Density Dist. 
for 28 Diam. Gap 

The corresponding results for a gap of 
14 diams. have been computed, and found to lie 
between those for 7 and 28 diameters. The 
lining up of particles, as observed in Figure 
1, is an extreme example of the "layering" 
reported by S&D. These investigators also 
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report that with a 7 diam. gap, the interior 
flow stresses agree reasonably well with 
extended kinetic theories. They do not give, 
however, results for greater gaps, where we 
find a trend away from the smaller gap behav- 
ior. It is pertinent to observe at this junc- 
ture that computer simulations of an infinite 
flow field do, by the very manner in which 
they mimic infinity in the direction trans- 
verse to the flow, tend to deny the possibili- 
ty of flow coagulat ion. 

The reader will have observed that the 
higher-density velocity data do not necessari- 
ly support an antisymmetric profile. Whereas 
at low densities, when the same gross boundary 
conditions are applied, all initial particle 
distributions appear to lead to the same mean 
flow, the same cannot be uneqivocally stated 
for higher densities. Lack of antisymmetry and 
of reproducibility has been observed. Further 
study should be made of this matter. S&D also 
report similar findings. 



Figure 11 Dependence of Normal Stress on 
Solids % and Gap 

The dependence of the normal stress on 
gap and solids content is shown in Figure 11. 
At low densities this stress varies nearly 
inversely to the gap. On the other hand, at 
the upper density limit, the stress is gener- 
ated almost totally by compression, and is 
independent of gap size. The corresponding 
values of the overall friction coefficient, 
shear/normal stress, are shown in Figure 12. 
There is a change in behavior with gap size. 
The curve for the largest gap exhibits the 
minimum observed for infinite media. 

To determine the influence of the inter- 
nal friction coefficient, CF, runs were made 
with a gap of 7 diams. Figure 13 shows the 
significant, but not dramatic, difference 
between the results for CF = 0 and CF = 0.3. 

w SLIDER BEARING ”: 

The second case treated numerically was 
an inclined "slider bearing” [For a discussion 
of this case, see ref. 5]. A uniform layer of 
granules is drawn by a moving plate into 
the wedge formed by that plate and another, of 
finite length, inclined to it. Figure 14 shows 





Figure 12 Overall Friction Coefficient as 
Function of Solids % and Gap 


the "slider” with an inlet gap of 20 granule 
diameters, an exit gap of 10 diams., and a 
total length of 20 diams. In every case, the 
thickness of the oncoming layer is such that 
it just fits into the exit gap, H2 (thus a 
thickness of 4 for an exit H2 = 5) . 



Figure 13 Dependence of Overall Friction 
Coeff.on Friction Coeff. between Granules 



Figure 14 Entrance Region of "Slider” 
with HI = 10 Diam. & H2 = 5 Diam. 

Figure 14 shows the entrance region. We 
see the particles rejected backwards as the 
slider plows its way over the oncoming carpet 
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slider plows its way over the oncoming carpet 
of grannies. In Fig. 15 we see granules burst 
out from the compression experienced in the 
tight gap. The velocity profiles at entrance, 
midway and exit are shown in Figure 16. They 
resemble those to be found in a classical, 
liquid-lubricated flat slider bearing. 



Figure 15 Exit Region of Slider 
with HI = 10 Diam. fc H2 = 5 Diam. 


The force distribution exerted by the 
granules upon the slider was not calculated, 


but the relative 

total longitudinal and verti 

cal forces 

were 

determined, as 

given below. 

Configuration 



HI 

H2 

D 

FX 

FY 

10 

5 

10 

2.58 

4.00 

10 

5 

20 

3.94 

8.11 

8 

8 

10 

2.36 

4.46 


The significant result here is that, within 
the range of parameters investigated, the load 
capacity does not depend on the slider incli- 
nation . Thus, although the velocity profiles 
may resemble those of classical lubrication, 
the load-carrying capacity does not. 



Figure 16 Velocity Dist. within "Slider" 
with HI = 10 Diam. & H2 = 5 Diam 
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28 particle diameters are investigated, with solid fractions ranging from 0.2 to 0.9. Among features observed are: slip 
flow at boundaries, coagulation at high densities, and gross fluctuations in surface stress. A videotape has been 
prepared to demonstrate the foregoing effects. 
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